﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "netn")]
    public static unsafe int netn(int n, double eps, double h, IntPtr x_ptr, IntPtr f_xa_ya_n_ptr, int inter)
    {
        double* x = (double*)x_ptr.ToPointer();
        f_xa_ya_n = Marshal.GetDelegateForFunctionPointer<delegatefunc_xa_ya_n>(f_xa_ya_n_ptr);

        return netn(n, eps, h, x, inter);
    }

    /// <summary>
    /// 非线性方程组拟牛顿法
    /// f(double[] x, double[] y, int n) 计算方程组各方程左端函数值的函数名。
    /// </summary>
    /// <param name="n">非方程组阶数。</param>
    /// <param name="eps">控制精度要求。</param>
    /// <param name="h">增量初值,h>0。返回时将被破坏。</param>
    /// <param name="x">x[n]存放初值。返回方程组实数解。</param>
    /// <param name="inter">最多迭代次数</param>
    /// <returns>函数返回实际迭代次数。若<0则表示因AZ=B奇异或beta=0而失败。本函数最大迭代次数为1000。</returns>
    public static unsafe int netn(int n, double eps, double h, double* x, int inter = 1000)
    {
        //控制h大小
        int i, j, k;
        double t, ep, z, beta, d;
        double* y = stackalloc double[n];
        double* a = stackalloc double[n * n];
        double* b = stackalloc double[n];

        k = 0;
        t = 0.1;
        ep = 1.0 + eps;
        while (ep >= eps)
        {
            f_xa_ya_n(x, b, n);
            ep = 0.0;
            for (i = 0; i <= n - 1; i++)
            {
                z = Math.Abs(b[i]);
                if (z > ep) ep = z;
            }
            if (ep >= eps)
            {
                k = k + 1;
                // 达到最大迭代次数未收敛
                if (k == inter)
                {
                    return (k);
                }
                for (j = 0; j <= n - 1; j++)
                {
                    z = x[j];
                    x[j] = x[j] + h;
                    f_xa_ya_n(x, y, n);
                    for (i = 0; i <= n - 1; i++)
                    {
                        a[i * n + j] = y[i];
                    }
                    x[j] = z;
                }
                // AZ=B奇异
                if (gauss(a, b, n) == 0)
                {
                    return (-1);
                }
                beta = 1.0;
                for (i = 0; i <= n - 1; i++)
                {
                    beta = beta - b[i];
                }

                // beta=0
                //if (Math.Abs(beta) + 1.0 == 1.0)
                if (EQ0(beta))
                {
                    return (-1);
                }
                d = h / beta;
                for (i = 0; i <= n - 1; i++)
                {
                    x[i] = x[i] - d * b[i];
                }
                h = t * h;
            }
        }
        return (k);
    }

    /*
      int main()
      { 
          int i, k;
          void netnf(double [],double [],int);
          double eps, h;
          double x[3]={1.0,1.0,1.0};
          h=0.1; eps=0.0000001; 
          k = netn(3,eps,h,x,netnf);
          cout <<"迭代次数 = " <<k <<endl;
          for (i=0; i<=2; i++)  cout <<"x(" <<i <<") = " <<x[i] <<endl;
          return 0;
      }
    // 方程组
      void netnf(double x[], double y[], int n)
      { y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-1.0;
        y[1]=2.0*x[0]*x[0]+x[1]*x[1]-4.0*x[2];
        y[2]=3.0*x[0]*x[0]-4.0*x[1]+x[2]*x[2];
        n=n;
        return;
      }
    */
}

